## Wild orangutans maintain sleep homeostasis through napping, counterbalancing socio-ecological 
    # factors that interfere with their sleep
## Analysis Script

# 0) Packages ----

library(lme4)
library(lmerTest)
library(multcomp)
library(performance)
library(effects)
library(ggeffects)
library(DHARMa)
library(colorspace)
library(lattice)
library(glmmTMB)
library(plotrix)

# 1) NIGHTTIME SLEEP PERIOD -----------------------------------------------------------------------------

# >> Load data ----

#setwd("")
df_night <- read.csv("NighttimeSleepPeriodDuration_data_fin.csv", sep = ";")
df_night$focalID <- as.factor(df_night$focalID)

# >> Checks ----

str(df_night)
summary(df_night)
head(df_night)

df_night$sunrise <- as.POSIXct(df_night$sunrise, format = "%d.%m.%Y %H:%M", tz = "UTC")
df_night$sunset <- as.POSIXct(df_night$sunset, format = "%d.%m.%Y %H:%M", tz = "UTC")

hist(df_night$FAI_Fruit)
hist(df_night$trj_DJL)
hist(df_night$CountPMAtRS_AllDist_InfCorr) 
hist(df_night$SleepDurationNextNight) 

# > Summary stats --------------------------------------------------------------------------------------

#sunrise and sunset
#take all of the times and paste them together with the same day (today, but any day is fine)
today_times_sunrise <- as.POSIXct(paste("2025-03-06", substr(df_night$sunrise, 12, 19)), tz = "UTC")
today_times_sunset <- as.POSIXct(paste("2025-03-06", substr(df_night$sunset, 12, 19)), tz = "UTC")

#calculate mean sunset and sunrise times
mean(today_times_sunrise) + (7*60*60) #add seven hours (to get to WIB tz)
sd(today_times_sunrise)/60 #divide by 60s to get minutes
min(today_times_sunrise) + (7*60*60) #add seven hours (to get to WIB tz)
max(today_times_sunrise) + (7*60*60) #add seven hours (to get to WIB tz)

mean(today_times_sunset)  + (7*60*60) #add seven hours (to get to WIB tz)
sd(today_times_sunset)/60 #divide by 60s to get minutes
min(today_times_sunset) + (7*60*60) #add seven hours (to get to WIB tz)
max(today_times_sunset) + (7*60*60) #add seven hours (to get to WIB tz)

# night nest start
mean(df_night$BedTimeDayBefore_Dec)
(mean(df_night$BedTimeDayBefore_Dec)-17)*60
min(df_night$BedTimeDayBefore_Dec)
(min(df_night$BedTimeDayBefore_Dec)-13)*60
max(df_night$BedTimeDayBefore_Dec)
(max(df_night$BedTimeDayBefore_Dec)-19)*60

# night nest end
mean(df_night$WakeUpTimeDayAfterDec)
(mean(df_night$WakeUpTimeDayAfterDec)-6)*60
min(df_night$WakeUpTimeDayAfterDec)
(min(df_night$WakeUpTimeDayAfterDec)-5)*60
max(df_night$WakeUpTimeDayAfterDec)
(max(df_night$WakeUpTimeDayAfterDec)-8)*60

# nighttime sleep period
mean(df_night$SleepDurationNextNight)
(mean(df_night$SleepDurationNextNight)-12)*60

sd(df_night$SleepDurationNextNight)
sd(df_night$SleepDurationNextNight)*60


# > Model ------------------------------------------------------------------------------------------------

glmmSleepDuration_Full<-lmer(SleepDurationNextNight  ~
                               +scale(trj_DJL)
                             +scale(FAI_Fruit)
                             +scale(CountPMAtRS_AllDist_InfCorr)
                             +scale(SleepDurationNightBefore)
                             +scale(rainfallMorningAfter)
                             #+scale(temperatureMaxNEXTMorning)
                             +scale(temperatureMinNEXTMorning)
                             +scale(TotalNapTime)
                             +scale(TBuildNightNest)
                             #+ComfortElement
                             +scale(NightDuration_Dec)
                             +scale(KcalIntake)
                             +MoonFraction
                             +Class
                             +(1|focalID),data=df_night)

glmmSleepDuration_Null<-lmer(SleepDurationNextNight  ~
                               +(1|focalID),data=df_night)

anova(glmmSleepDuration_Full,glmmSleepDuration_Null)

summary(glmmSleepDuration_Full) 
cftest(glmmSleepDuration_Full) 
conf_intervals_wald_DJL_nightsleep <- confint(glmmSleepDuration_Full, method = "Wald")
print(conf_intervals_wald_DJL_nightsleep)

summary(glht(glmmSleepDuration_Full, linfct = mcp(Class = "Tukey")), test = adjusted("holm"))
check_collinearity(glmmSleepDuration_Full)
r2_nakagawa(glmmSleepDuration_Full)
qqmath(glmmSleepDuration_Full)
plot(glmmSleepDuration_Full)
plot(df_night$LengthActivitySameDay_Dec~df_night$trj_DJL)

# >> Effect sizes ----------------------------------------------------------------------------------------


#assoc
summary(df_night$CountPMAtRS_AllDist_InfCorr)
Preds<-ggeffect(glmmSleepDuration_Full, terms = "CountPMAtRS_AllDist_InfCorr[0:5]")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change) *60 # in mins
#sd(Preds_Change)

#djl
dists <- seq(100,2400,100)
summary(df_night$trj_DJL)
Preds<-ggeffect(glmmSleepDuration_Full, terms = "trj_DJL[dists]")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change) *60 #in mins per 100m


#night before sleep
summary(df_night$SleepDurationNightBefore)
Preds<-ggeffect(glmmSleepDuration_Full, terms = "SleepDurationNightBefore[9:17]")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change)  *60 #in mins per hour


#min temp
summary(df_night$temperatureMinNEXTMorning)
Preds<-ggeffect(glmmSleepDuration_Full, terms = "temperatureMinNEXTMorning[20:26]")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change)  *60 #in mins per hour


# > Plots ----------------------------------------------------------------------------------------


#Preps for plots

ColAdFem<-"#995D81"
ColFlMale<-"#006BA6"
ColMoth<-"#A30000"
ColUFM<-"#70A288"

ColAdFem<-adjust_transparency(ColAdFem, alpha = 0.5)
ColUFM<-adjust_transparency(ColUFM, alpha = 0.5)
ColFlMale<-adjust_transparency(ColFlMale, alpha = 0.5)
ColMoth<-adjust_transparency(ColMoth, alpha = 0.5)

Symb2AdFem<-15
Symb2FlMale<-23
Symb2Moth<-16
Symb2UFM<-17

Color<-data.frame(cbind(1:4,
                        c(ColAdFem,ColFlMale,ColMoth,ColUFM),
                        c(Symb2AdFem,Symb2FlMale,Symb2Moth,Symb2UFM),
                        c("Adult Female", "Flanged Male", "Mother", "Unflanged Male")))

colnames(Color)<-c("ClassNr","ColorClass","Pch2Class", "Class")
dataModelSleepDur<- merge(df_night,Color,by="Class")

str(dataModelSleepDur)
head(dataModelSleepDur)

#Figure2 

#prepping plot layout
par(mar=c(5.1, 6.1, 3.1, 2.1))

layout(matrix(c(1,2,
                1,2,
                1,2,
                1,2,
                3,4,
                3,4,
                3,4,
                3,4,
                5,5), nrow=9, byrow=TRUE))
layout.show(n=5)

#panel a
#par(mfrow=c(1,1))

plot(SleepDurationNextNight~jitter(CountPMAtRS_AllDist_InfCorr,1.5),
     #xlim=c(0.4,7.1),
     data=dataModelSleepDur,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     axes=F,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="")

title(ylab="Night sleep period duration (h)", line=4, cex.lab=2)
title(xlab="No. of nighttime association partners", line=4, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=10,to=17,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=-1,to=12,by=1)))

PredsCountPMAtRS_AllDist_InfCorr<-ggeffect(glmmSleepDuration_Full, terms = "CountPMAtRS_AllDist_InfCorr")
lines(PredsCountPMAtRS_AllDist_InfCorr$predicted~PredsCountPMAtRS_AllDist_InfCorr$x, lty=1, lwd=3,col="grey30")
str(PredsCountPMAtRS_AllDist_InfCorr)


lower_ci_1a <- PredsCountPMAtRS_AllDist_InfCorr$conf.low
upper_ci_1a <- PredsCountPMAtRS_AllDist_InfCorr$conf.high


ylow_1a<-PredsCountPMAtRS_AllDist_InfCorr[,"conf.low"]
yhigh_1a<-rev(PredsCountPMAtRS_AllDist_InfCorr[,"conf.high"])

Dur1a <- PredsCountPMAtRS_AllDist_InfCorr[,1]

xcoord1a<-c((Dur1a),rev(Dur1a))
ycoord1a<-c(ylow_1a,yhigh_1a)

polygon(x=xcoord1a,
        y=ycoord1a, 
        col = "#80808040", border = NA)

corner.label(label='A',figcorner=T,cex=2)

#panel b
#par(mfrow=c(1,1))
par(family="sans")

plot(SleepDurationNextNight~jitter(trj_DJL,2),
     data=dataModelSleepDur,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     axes=F,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="")

title(ylab="Night sleep period duration (h)", line=4, cex.lab=2)
title(xlab="Day journey length (m)", line=4, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=10,to=17,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=-200,to=2500,by=200)))

PredsDJL<-ggeffect(glmmSleepDuration_Full, terms = "trj_DJL")
lines(PredsDJL$predicted~PredsDJL$x, lty=1, lwd=3,col="grey30")
str(PredsDJL)

lower_ci_1b <- PredsDJL$conf.low
upper_ci_1b <- PredsDJL$conf.high


ylow_1b<-PredsDJL[,"conf.low"]
yhigh_1b<-rev(PredsDJL[,"conf.high"])

Dur1b <- PredsDJL[,1]

xcoord1b<-c((Dur1b),rev(Dur1b))
ycoord1b<-c(ylow_1b,yhigh_1b)

polygon(x=xcoord1b,
        y=ycoord1b, 
        col = "#80808040", border = NA)


corner.label(label='B',figcorner=T,cex=2)

#panel c
#par(mfrow=c(1,1))

plot(SleepDurationNextNight~jitter(SleepDurationNightBefore,2),
     data=dataModelSleepDur,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     axes=F,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="")

title(ylab="Night sleep period duration (h)", line=4, cex.lab=2)
title(xlab="Preceeding night sleep period duration (h)", line=4, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=10,to=17,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=9,to=17,by=1)))

PredsSleepDurationNightBefore<-ggeffect(glmmSleepDuration_Full, terms = "SleepDurationNightBefore")
lines(PredsSleepDurationNightBefore$predicted~PredsSleepDurationNightBefore$x, lty=1, lwd=3,col="grey30")

str(PredsSleepDurationNightBefore)

lower_ci_1c <- PredsSleepDurationNightBefore$conf.low
upper_ci_1c <- PredsSleepDurationNightBefore$conf.high


ylow_1c<-PredsSleepDurationNightBefore[,"conf.low"]
yhigh_1c<-rev(PredsSleepDurationNightBefore[,"conf.high"])

Dur1c <- PredsSleepDurationNightBefore[,1]

xcoord1c<-c((Dur1c),rev(Dur1c))
ycoord1c<-c(ylow_1c,yhigh_1c)

polygon(x=xcoord1c,
        y=ycoord1c, 
        col = "#80808040", border = NA)


corner.label(label='C',figcorner=T,cex=2)

#panel d
#par(mfrow=c(1,1))

plot(SleepDurationNextNight~temperatureMinNEXTMorning,
     data=dataModelSleepDur,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     axes=F,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="")

title(ylab="Night sleep period duration (h)", line=4, cex.lab=2)
title(xlab="Minimum nighttime temperature (°C)", line=4, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=10,to=17,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=19,to=27,by=1)))

PredsMinTemp<-ggeffect(glmmSleepDuration_Full, terms = "temperatureMinNEXTMorning")
lines(PredsMinTemp$predicted~PredsMinTemp$x, lty=1, lwd=3,col="grey30")
str(PredsMinTemp)

lower_ci_1d <- PredsMinTemp$conf.low
upper_ci_1d <- PredsMinTemp$conf.high


ylow_1d<-PredsMinTemp[,"conf.low"]
yhigh_1d<-rev(PredsMinTemp[,"conf.high"])

Dur1d <- PredsMinTemp[,1]

xcoord1d<-c((Dur1d),rev(Dur1d))
ycoord1d<-c(ylow_1d,yhigh_1d)

polygon(x=xcoord1d,
        y=ycoord1d, 
        col = "#80808040", border = NA)

corner.label(label='D',figcorner=T,cex=2)


#adding legend
#par(mfrow=c(1,1))

par(mar=c(0,0,0,0))
plot(NULL ,xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("center",legend=c("Adult females","Flanged Males","Mothers","Unflanged Males"),
       col=c(ColAdFem,ColFlMale,ColMoth,ColUFM),
       pch=c(Symb2AdFem,18,Symb2Moth,Symb2UFM),
       pt.lwd=1,
       cex=2,
       pt.cex=c(1.5,1.9,1.6,1.6),horiz=TRUE)





# 2a) DAYTIME NAP PERIODS --------------------------------------------------------------------------------------------

# >> Load data ----

df_day <- read.csv("DaytimeCumulativeNapDuration_data_fin.csv", sep = ";")
df_day$focalID <- as.factor(df_day$focalID)

# >> Checks ----

str(df_day)
head(df_day)

df_day$Class <- as.factor(df_day$Class)


# > Summary stats -------------------------------------------------------------------------------------

# days without any naps
sum(df_day$TotalNapTime == 0)
268/455
455-268
187/455*100 #percentage of follows with at least 1 daynest


#cumulative duration of naps
mean(df_day$TotalNapTime[which(df_day$TotalNapTime >0)])
sd(df_day$TotalNapTime[which(df_day$TotalNapTime >0)])


# > Model ---------------------------------------------------------------------------------------------
glmmDayNapDuration_Full_PoisZeroInfl2<-glmmTMB(TotalNapTime ~
                                                 scale(trj_DJL)
                                               +scale(FAI_Fruit)
                                               +scale(AverageAssPartners) 
                                               +scale(SleepDurationNightBefore) 
                                               +scale(MaxTempDuringDay)
                                               +scale(KcalIntake)
                                               +RainfallDuringDayYN
                                               +Class
                                               +scale(KcalIntake)
                                               +(1|focalID),
                                               zi=~log(LentghtActivitySameDay_Dec), 
                                               family = poisson,
                                               data=df_day)



glmmDayNapDuration_Null_PoisZeroInfl2<-glmmTMB(TotalNapTime ~
                                                 +(1|focalID),
                                               zi=~log(LentghtActivitySameDay_Dec), 
                                               family = poisson,
                                               data=df_day)

anova(glmmDayNapDuration_Full_PoisZeroInfl2,glmmDayNapDuration_Null_PoisZeroInfl2) # full model sig better than null model
summary(glmmDayNapDuration_Full_PoisZeroInfl2)
conf_intervals_wald_DJL <- confint(glmmDayNapDuration_Full_PoisZeroInfl2, method = "Wald")
print(conf_intervals_wald_DJL)
testDispersion(glmmDayNapDuration_Full_PoisZeroInfl2)
testZeroInflation(glmmDayNapDuration_Full_PoisZeroInfl2)
summary(glht(glmmDayNapDuration_Full_PoisZeroInfl2, linfct = mcp(Class = "Tukey")), test = adjusted("holm"))
r2_nakagawa(glmmDayNapDuration_Full_PoisZeroInfl2)
check_collinearity(glmmDayNapDuration_Full_PoisZeroInfl2)
residglmmDayNapDuration_Full_PoisZeroInfl2 <- simulateResiduals(glmmDayNapDuration_Full_PoisZeroInfl2)
plot(residglmmDayNapDuration_Full_PoisZeroInfl2)


# >> Effect sizes --------------------------------------------------------------------------------------------


#sleep duration night before
summary(df_day$SleepDurationNightBefore)
Preds_sleepdurnightbefore<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "SleepDurationNightBefore[9:17]")
Preds_sleepdurnightbefore_Decrease <- Preds_sleepdurnightbefore$predicted[2:9] - Preds_sleepdurnightbefore$predicted[1:8]
mean(Preds_sleepdurnightbefore_Decrease)
sd(Preds_sleepdurnightbefore_Decrease)
#beta for scale(sleepdurationnightbefore) is -0.095431
Preds_sleepdurnightbefore$predicted[1:8]/Preds_sleepdurnightbefore$predicted[2:9]
1-(Preds_sleepdurnightbefore$predicted[2:9]/Preds_sleepdurnightbefore$predicted[1:8])


# social partners
summary(df_day$AverageAssPartners)
Preds<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "AverageAssPartners[0:7]")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change)
sd(Preds_Change)
1- (Preds$predicted[2:8]/Preds$predicted[1:7])


# temp
summary(df_day$MaxTempDuringDay)
Preds<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "MaxTempDuringDay[23:38]")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change)
sd(Preds_Change)
1- (Preds$predicted[2:16]/Preds$predicted[1:15])

# calories
cal_seq <- seq(from = 730, to = 14556, by = 500)
summary(df_day$KcalIntake)
Preds<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "KcalIntake[cal_seq]")
Preds$predicted[1:27]/Preds$predicted[2:28]
1- (Preds$predicted[2:28]/Preds$predicted[1:27])

# rain
summary(df_day$RainfallDuringDayYN)
Preds<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "RainfallDuringDayYN")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
mean(Preds_Change)
sd(Preds_Change)

# AS class
summary(df_day$Class)
Preds<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "Class")
Preds_Change <- Preds$predicted[2:nrow(Preds)] - Preds$predicted[1:(nrow(Preds)-1)]
Preds

#djl
dists <- seq(100,2400,100)
summary(df_day$trj_DJL)
Preds<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "trj_DJL[dists]")
Preds$predicted[1:23]/Preds$predicted[2:24]
1- (Preds$predicted[2:24]/Preds$predicted[1:23])

# 2b) NUMBER OF NAPS ---------------------------------------------------------------------------------

# >> Load data ----

df_daynests <- read.csv("AverageDurationOfNapsPerDay_data_fin.csv", sep = ";")

# >> Checks ----

summary(df_daynests)
df_daynests$focalID<-as.factor(df_daynests$focalID)
str(df_daynests)
head(df_daynests)

# > Summary stats ------------------------------------------------------------------------------------

mean(df_daynests$NrDayNests)
sd(df_daynests$NrDayNests)

# > Model ---------------------------------------------------------------------------------------------

glmmNapHomeo_Full<-lmer(log(AverageNapTimePerDayNest+1) ~ NrDayNests 
                        +(1|focalID),data=df_daynests)

glmmNapHomeo_Null<-lmer(AverageNapTimePerDayNest ~
                          +(1|focalID),data=df_daynests)

anova(glmmNapHomeo_Full,glmmNapHomeo_Null)

summary(glmmNapHomeo_Full)
cftest(glmmNapHomeo_Full)

conf_intervals_wald_DJL_naphomeo <- confint(glmmNapHomeo_Full, method = "Wald")
print(conf_intervals_wald_DJL_naphomeo)


plot(glmmNapHomeo_Full)
qqnorm(resid(glmmNapHomeo_Full))
qqline(resid(glmmNapHomeo_Full))
r2_nakagawa(glmmNapHomeo_Full)


# > Plots ----------------------------------------------------------------------------------------------


#Preps: adding color and symbol data to the DF

ColAdFem<-"#995D81"
ColFlMale<-"#006BA6"
ColMoth<-"#A30000"
ColUFM<-"#70A288"

ColAdFem<-adjust_transparency(ColAdFem, alpha = 0.5)
ColUFM<-adjust_transparency(ColUFM, alpha = 0.5)
ColFlMale<-adjust_transparency(ColFlMale, alpha = 0.5)
ColMoth<-adjust_transparency(ColMoth, alpha = 0.5)


Symb2AdFem<-15
Symb2FlMale<-23
Symb2Moth<-16
Symb2UFM<-17

Color<-data.frame(cbind(1:4,
                        c(ColAdFem,ColFlMale,ColMoth,ColUFM),
                        c(Symb2AdFem,Symb2FlMale,Symb2Moth,Symb2UFM),
                        c("Adult Female", "Flanged Male", "Mother", "Unflanged Male")))

colnames(Color)<-c("ClassNr","ColorClass","Pch2Class", "Class")

df_day <- merge(df_day,Color,by="Class")
head(df_day)
summary(df_day)

df_daynests <- merge(df_daynests,Color,by="Class")
summary(df_daynests)
head(df_day,100)

#Figure 2

#prepping plot layout
par(mar=c(5.1, 6.1, 3.1, 2.1))
layout(matrix(c(1,2,3,
                1,2,3,
                1,2,3,
                1,2,3,
                4,5,6,
                4,5,6,
                4,5,6,
                4,5,6,
                7,7,7), nrow=9, byrow=TRUE))
layout.show(n=7)


#panel a
#par(mfrow=c(1,1))
par(family="sans")

plot(TotalNapTime/60~SleepDurationNightBefore,
     data=df_day,
     ylim=c(0,5.5),
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="",
     axes=F)

title(ylab="Cumulative daily nap period duration (h)", line=4, cex.lab=2)
title(xlab="Preceeding night sleep period duration (h)", line=3, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=-1,to=6,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=9,to=17,by=1)))

PredsDJL_TotalNap<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "SleepDurationNightBefore")
lines((PredsDJL_TotalNap$predicted)/60~ PredsDJL_TotalNap$x, lty=1, lwd=3,col="grey30")

lower_ci_2a <- PredsDJL_TotalNap$conf.low
upper_ci_2a <- PredsDJL_TotalNap$conf.high
str(PredsDJL_TotalNap)

ylow_2a<-PredsDJL_TotalNap[,"conf.low"]/60
yhigh_2a<-rev(PredsDJL_TotalNap[,"conf.high"])/60

Dur <- PredsDJL_TotalNap[,1]

xcoord<-c((Dur),rev(Dur))
ycoord<-c(ylow_2a,yhigh_2a)

polygon(x=xcoord,
        y=ycoord, 
        col = "#80808040", border = NA)

corner.label(label="A",figcorner=T,cex=2)


#panel b
#par(mfrow=c(1,1))
plot(TotalNapTime/60~AverageAssPartners,
     ylim=c(0,5.5),
     data=df_day,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="",
     axes=F)

title(ylab="Cumulative daily nap period duration (h)", line=3.8, cex.lab=2)
title(xlab="Average no. of daytime association partners", line=3, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=-1,to=6,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=-1,to=7,by=1)))

PredsAverageAssPartners_TotalNap<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "AverageAssPartners")
lines(PredsAverageAssPartners_TotalNap$predicted/60~ PredsAverageAssPartners_TotalNap$x, 
      lty=1, lwd=3,col="grey30")
str(PredsAverageAssPartners_TotalNap)

lower_ci_2b <- PredsAverageAssPartners_TotalNap$conf.low
upper_ci_2b <- PredsAverageAssPartners_TotalNap$conf.high

ylow_2b<-PredsAverageAssPartners_TotalNap[,"conf.low"]/60
yhigh_2b<-rev(PredsAverageAssPartners_TotalNap[,"conf.high"])/60

Dur2b <- PredsAverageAssPartners_TotalNap[,1]

xcoord_2b<-c((Dur2b),rev(Dur2b))
ycoord_2b<-c(ylow_2b,yhigh_2b)

polygon(x=xcoord_2b,
        y=ycoord_2b, 
        col = "#80808040", border = NA)



corner.label(label='B',figcorner=T,cex=2)

#panel c
#par(mfrow=c(1,1))
plot(TotalNapTime/60~KcalIntake,
     ylim=c(0,5.5),
     data=df_day,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     pch=as.numeric(Pch2Class),
     axes=F,
     ylab="",
     xlab="")

title(ylab="Cumulative daily nap period duration (h)", line=3.8, cex.lab=2)
title(xlab="Energy intake (Kcal)", line=3, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=-1,to=6,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=0,to=16000,by=2000)))

PredsKcalIntake_TotalNap<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "KcalIntake")
lines(PredsKcalIntake_TotalNap$predicted/60~ PredsKcalIntake_TotalNap$x, 
      lty=1, lwd=3,col="grey30")

str(PredsKcalIntake_TotalNap)

lower_ci_2c <- PredsKcalIntake_TotalNap$conf.low
upper_ci_2c <- PredsKcalIntake_TotalNap$conf.high

ylow_2c<-PredsKcalIntake_TotalNap[,"conf.low"]/60
yhigh_2c<-rev(PredsKcalIntake_TotalNap[,"conf.high"])/60

Dur2c <- PredsKcalIntake_TotalNap[,1]

xcoord_2c<-c((Dur2c),rev(Dur2c))
ycoord_2c<-c(ylow_2c,yhigh_2c)

polygon(x=xcoord_2c,
        y=ycoord_2c, 
        col = "#80808040", border = NA)


corner.label(label='C',figcorner=T,cex=2)

#panel d
#par(mfrow=c(1,1))
boxplot(TotalNapTime/60~RainfallDuringDayYN,
        ylim=c(0,5.5),
        data=df_day,
        cex=1.5,
        cex.lab=2,
        cex.axis=1.8,
        lwd=2,
        yaxt = "n",
        ylab="",
        xlab="Rainfall during the day",
        outline=F,
        cex.names=0.2,
        border="white",
        names=c("No","Yes"),
        col="white")
axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=-1,to=5.5,by=1)))
title(ylab="Cumulative daily nap period duration (h)", line=4, cex.lab=2)


df_dayRain<-df_day[df_day$RainfallDuringDayYN==1,]
df_dayNoRain<-df_day[df_day$RainfallDuringDayYN==0,]

points(jitter(rep(1,length(df_dayNoRain$TotalNapTime)),factor=10),
       df_dayNoRain$TotalNapTime/60, lwd=2,
       col=df_dayNoRain$ColorClass,bg=ColFlMale,
       pch=as.numeric(df_dayNoRain$Pch2Class))
points(jitter(rep(2,length(df_dayRain$TotalNapTime)),factor=5),
       df_dayRain$TotalNapTime/60, lwd=2,
       col=df_dayRain$ColorClass,bg=ColFlMale,
       pch=as.numeric(df_dayRain$Pch2Class))

#adding significance indicator
segments(x0 = 1, y0 = 5.3,
         x1 = 2, y1 = 5.3,
         lwd=2,col="black")
text(1.5,5.45,"***",cex=2,col="black")

#adding model predictions
PredsRain_TotalNap<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, 
                             terms = "RainfallDuringDayYN")

segments(x0 = 0.75, y0 = PredsRain_TotalNap[1,2]/60,
         x1 = 1.25, y1 = PredsRain_TotalNap[1,2]/60,
         lwd=3,col="grey30",lty=1)

segments(x0 = 0.75, y0 = PredsRain_TotalNap[1,"conf.low"]/60,
         x1 = 1.25, y1 = PredsRain_TotalNap[1,"conf.low"]/60,
         lwd=2,col="grey30")

segments(x0 = 0.75, y0 = PredsRain_TotalNap[1,"conf.high"]/60,
         x1 = 1.25, y1 = PredsRain_TotalNap[1,"conf.high"]/60,
         lwd=2, col="grey30")

segments(x0 = 1, y0 = PredsRain_TotalNap[1,"conf.low"]/60,
         x1 = 1, y1 = PredsRain_TotalNap[1,"conf.high"]/60,
         lwd=2, col="grey30")

segments(x0 = 1.75, y0 = PredsRain_TotalNap[2,2]/60,
         x1 = 2.25, y1 = PredsRain_TotalNap[2,2]/60,
         lwd=3,col="grey30",lty=1)

segments(x0 = 1.75, y0 = PredsRain_TotalNap[2,"conf.low"]/60,
         x1 = 2.25, y1 = PredsRain_TotalNap[2,"conf.low"]/60,
         lwd=2,col="grey30")

segments(x0 = 1.75, y0 = PredsRain_TotalNap[2,"conf.high"]/60,
         x1 = 2.25, y1 = PredsRain_TotalNap[2,"conf.high"]/60,
         lwd=2,col="grey30")

segments(x0 = 2, y0 = PredsRain_TotalNap[2,"conf.low"]/60,
         x1 = 2, y1 = PredsRain_TotalNap[2,"conf.high"]/60,
         lwd=2,col="grey30")

corner.label(label='D',figcorner=T,cex=2)

#panel e
#par(mfrow=c(1,1))

plot(TotalNapTime/60~MaxTempDuringDay,
     ylim=c(0,5.5),
     data=df_day,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     pch=as.numeric(Pch2Class),
     axes=F,
     ylab="",
     xlab="")

title(ylab="Cumulative daily nap period duration (h)", line=4, cex.lab=2)
title(xlab="Maximum daytime temperature (°C)", line=4, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=-1,to=6,by=1)))
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=15,to=40,by=2)))

PredsTemp_TotalNap<-ggeffect(glmmDayNapDuration_Full_PoisZeroInfl2, terms = "MaxTempDuringDay")
lines(PredsTemp_TotalNap$predicted/60~ PredsTemp_TotalNap$x,
      lty=1, lwd=3,col="grey30")

str(PredsTemp_TotalNap)

lower_ci_2e <- PredsTemp_TotalNap$conf.low
upper_ci_2e <- PredsTemp_TotalNap$conf.high

ylow_2e<-PredsTemp_TotalNap[,"conf.low"]/60
yhigh_2e<-rev(PredsTemp_TotalNap[,"conf.high"])/60

Dur2e <- PredsTemp_TotalNap[,1]

xcoord_2e<-c((Dur2e),rev(Dur2e))
ycoord_2e<-c(ylow_2e,yhigh_2e)

polygon(x=xcoord_2e,
        y=ycoord_2e, 
        col = "#80808040", border = NA)


corner.label(label='E',figcorner=T,cex=2)

#panel f
#par(mfrow=c(1,1))

plot(AverageNapTimePerDayNest~jitter(NrDayNests,1),
     xlim=c(0.5,4.5),
     ylim=c(-0.2,3.5),
     data=df_daynests,
     col=ColorClass,
     bg=ColFlMale,
     cex=1.5,
     cex.lab=2,
     cex.axis=2,
     lwd=2,
     pch=as.numeric(Pch2Class),
     ylab="",
     xlab="",
     axes=F)

title(ylab="Average nap period duration (h)", line=4, cex.lab=2)
title(xlab="Number of naps per day", line=4, cex.lab=2)

axis(side=2,las=2,mgp=c(3,1,0),cex.axis=1.8,at=c(seq(from=-1,to=6,by=1)),outer=F)
axis(side=1,las=1,mgp=c(3,1.2,0),cex.axis=1.8,at=c(seq(from=1,to=5,by=1),outer=F))

PredsNapHomeo<-ggeffect(glmmNapHomeo_Full, terms = "NrDayNests")
lines(PredsNapHomeo$predicted~PredsNapHomeo$x, lty=1, lwd=3,col="grey30")

str(PredsNapHomeo)

lower_ci_2f <- PredsNapHomeo$conf.low
upper_ci_2f <- PredsNapHomeo$conf.high

ylow_2f<-PredsNapHomeo[,"conf.low"]
yhigh_2f<-rev(PredsNapHomeo[,"conf.high"])

Dur2f <- PredsNapHomeo[,1]

xcoord_2f<-c((Dur2f),rev(Dur2f))
ycoord_2f<-c(ylow_2f,yhigh_2f)

polygon(x=xcoord_2f,
        y=ycoord_2f, 
        col = "#80808040", border = NA)


corner.label(label='F',figcorner=T,cex=2)


#adding legend
#par(mfrow=c(1,1))
par(mar=c(0,0,0,0))
plot(NULL ,xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("center",legend=c("Adult females","Flanged Males","Mothers","Unflanged Males"),
       col=c(ColAdFem,ColFlMale,ColMoth,ColUFM),
       pch=c(Symb2AdFem,18,Symb2Moth,Symb2UFM),
       pt.lwd=1,
       cex=2,
       pt.cex=c(1.5,1.9,1.6,1.6),horiz=TRUE)



# 3) FOLLOW UPS --------------------------------------------------------------------------------------


par(mfrow=c(1,1))

# > A) Compare association size at night nest vs during the day ----

str(df_night)
boxplot(df_night$AverageAssPartners,df_night$CountPMAtRS_AllDist_InfCorr)
t.test(df_night$AverageAssPartners,df_night$CountPMAtRS_AllDist_InfCorr)

mean(df_night$AverageAssPartners)
mean(df_night$CountPMAtRS_AllDist_InfCorr)

sd(df_night$AverageAssPartners)
sd(df_night$CountPMAtRS_AllDist_InfCorr)



# > B) Relationship between length of active period and DJL ----

glmmTravel_Full<-lmer(LengthActivitySameDay_Dec  ~
                        +scale(trj_DJL)
                      +Class
                      +(1|focalID),data=df_night)
summary(glmmTravel_Full) 
cftest(glmmTravel_Full)
conf_intervals_wald_DJL_Travel <- confint(glmmTravel_Full, method = "Wald")
print(conf_intervals_wald_DJL_Travel)
summary(glht(glmmTravel_Full, linfct = mcp(Class = "Tukey")), test = adjusted("holm"))
r2_nakagawa(glmmTravel_Full)

plot(LengthActivitySameDay_Dec  ~
       +scale(trj_DJL),data=df_night)

